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Abstract The Liege Oscillation code can be used as a stand- 
alone program or as a library of subroutines that the user 
calls from a Fortran main program of his own to compute 
radial and non-radial adiabatic oscillations of stellar models. 
We describe the variables and the equations used by the pro- 
gram and the methods used to solve them. A brief account is 
given of the use and the output of the program. 
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1 Introduction 

The Liege oscillation code (OSC) has been developed in 
the early 70s for computing adiabatic pulsations of speri- 
cally symmetric stars (no rotation nor magnetic field). It has 
gone through minor updates and is still presently in use in 
the asteroseismology group of the Liege Institute of Astro- 
physics and Geophysics. Besides the frequencies and eigen- 
functions, it produces also the coefficients needed to com- 
pute the first order rotational frequency splitting for a rigid 
rotation and the kernels needed to compute the splitting when 
a non rigid rotation is considered. 



2 Stellar models 

A stellar model is input to OSC as a table describing a few 
physical quantities at discrete points of the star, ordered from 
the centre to the surface. In theory, only two functions are 
necessary to compute stellar oscillations, for instance p (r) 
and T\(j), the density and the first adiabatic exponent in 
terms of the radius. However, for OSC, the model file must 
give at each point the values of the radius r, the mass m in 
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the spere of radius r, the total pressure P, the density p and 
the first adiabatic exponent T\ . It is clear that m and P could 
have been computed from the other quantities. The program 
does not require the Brunt- Vaisala frequency, often poorly 
computed by evolution codes. 

If the first point is not at the centre, OSC computes the 
oscillations of the given envelope with a rigid boundary con- 
dition at the bottom. Of course, when neglecting the oscil- 
latory behaviour of the core, great attention must be paid to 
the physical meaning of the output of the program. 

The outer boundary conditions will be applied at the last 
point, considered as the surface of the star. 

Inside the code, the model is described by the follow- 
ing five dimensionless quantities: x = r/R, q/x^ (with q = 
m/M), RP/GMp, AnR^p /M and Fy, where R and M denote 
the radius and the total mass of the star. 



3 Stellar oscillations 

3.1 Oscillation modes 

The small perturbations of a spherical star without rotation 
or magnetic field may be described as a superposition of nor- 
mal modes of oscillation which are the solutions of a linear 
boundary eigenvalue problem. These normal modes may be 
indexed by three integers k, I and m. Index k is loosely re- 
lated to the number of nodes of the radial displacement. In- 
dices I and m are the usual indices of the spherical function 
Y(,„{d,^) describing the angular dependence. Index £ may 
take any null or positive integer value and in may take 2£+l 
values between —i and +£. 

In the following description we use the notation dX and 
X' for the Lagrangian and Eulerian perturbation of any quan- 
tity X and a for the angular frequency. We often use the 
dimensionless angular frequency (0 = OTdyn, where the dy- 
namical time Tdyn is defined as \JWJgM. 
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3.2 Oscillation equations 



3.2.2 Non radial oscillations. 



The theory of stellar oscillation has been developed in a 
numb er of textbooks. We are the m ost familiar wit h the pa- 
per of iLedoux and WalravenI (Il958h and the books of lUnno eT 
(119791) and lCox. C198Q) . We will just write the needed equa- 



tions in the form they are implemented in our code. 
3.2.1 Radial oscillations. 

In the case of radial oscillations {£ = jn = 0), the equation of 
Poisson can be integrated and the perturbation of the gravi- 
tational potential eliminated. The differential system is then 
reduced to order two. We describe a normal mode with two 
functions Y{x) and Z{x). Disregarding an arbitrary phase, 
the displacement dr and the lagrangian perturbation of the 
pressure dP are written in terms of 7 (x) and Z{x) in the fol- 
lowing way. 



In the case of non radial oscillations {£ 0), the differential 
system is of order four. We describe a normal mode with four 
ailunctions Y{x), Z{x),U (x) and V{x). These functions as well 
as the frequency do not depend upon index m {2£ + 1 -fold 
degeneracy). The displacement reads 

a{r)YUe,(^)er + b{r){ '^\^!'^K e 



5r = 3i{a{r)e~""e,.} = V4n3i{a{r)Yooie,(^)e-""er} , (1) 

where e,- is a unit vector in the radial direction. Near the 
centre, a(r) oc r and may be written 

a{r)/r = Y{x) or a{r) / R = xY (x) . (2) 
In a similar way, 

5P/P = 3i{Zix)e-''"} = V47r9t{Z(jc)yoo(0,(/))e-'^'}. (3) 
Now, the differential equations read 



dx X Fix ' 

dZ / a ,\ GMp GMp q ^ 



dx 



RP 



RP x3 



(4) 

(5) 



These equations must be completed by the boundary condi- 
tions. At the centre, the regularity of the solution is ensured 
by the condition 



3riy+z = o at x = 0. 



(6) 



When an envelope model is given, the condition at the 
centre is replaced by a condition at the bottom of the enve- 
lope. 



y = 0. 



(7) 



At the surface, we generally apply a condition deduced 
from the vanishing of the pressure. In this case, the coeffi- 
cient GMp /RP in equation (|5]l tends to infinity and the reg- 
ularity of the solution requires that 



a^ + oAy + %rZ = o. 

x' I x^ 



(8) 



Note that if R is the value of r at the last point, x = qjy? — 1 
at the surface but this is not mandatory. The user can choose 
to apply the more usual condition 



5P = or Z = 0. 



(9) 



5r=\/47r9t 

, 1 ^yft„(0,</)) 



sin Q 



60 



(10) 



where e^, Ce and form the usual local cartesian basis of 
spherical coordinates. Near the centre, a{r) and b{r) oc 
and are written 

a(r)//?=/-'F(jc), 



b{r)/R=- 



J- 1 



(11) 

(12) 



The Lagrangian perturbation of pressure 5P and the Eule- 
rian perturbation of the gravitational potential 4'' are given 
by 

^/An'^{x^Z{x)Yi,„{e,(^)e-''"}, 



P 

Rd)' , 

— = V4^3i{/u{x)Yi„,{e,(^)e-'"}, 

UM 



Y(x) -Y(x) 



(13) 
(14) 

(15) 



GM dr \ 

xy|„(0,(/.)e-'^'|. 

Of course, the solution of a linear problem may be multi- 
plied by an arbitrary factor and the v^jr factor in the above 
expressions may be dropped. We have put it there for aes- 
thetic reasons because the spherical functions are normal- 
ized in such a way that 



4k 



\Yi„\^d^ = \. 



(16) 



With this V 471 factor, the time-average kinetic energy of a 
mode is given by 

Eun= J ^pv^dV = ^/ [a^+eie+l)b^]4nr^pdr 

= — [a^+£{£+\)b^]dm, (17) 

without any n factor. 

With these definitions, the oscillation equations read 

d^_l_ 
dx 



-Y 



X 
X 

dZ GMp ( 



RP 



GMp 



Z + U 



(18) 
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1) q 



RP 

GMp~ 



-Z + U 



dU 
dx 
dV_ 

dx 



X 



M 

(eu - V) 



■Y-IU 



xa)"- 



M 



4y 



RP 
GMp 



Z + U 



(19) 



(20) 



(21) 



These equations have been published byl Bourv et al.l(ll975h . 
but their equation (9) has been affected by a typo (an extra 
factor £). 

The regularity of the solution at the centre imposes two 
conditions at x = 0. 



£ 



x^ 



AnR^i. 
M 



RP 

" GMp 

-m. 



Z + U 



(22) 



(23) 



In the case of an envelope, the bottom of the envelope is 
supposed to behave as a rigid boundary. 



Y=0. 



(24) 



There are no movement nor density perturbations in the core, 
where the perturbation of the gravitational potential obeys a 
Laplace equation and has a simple analytical expression. We 
thus require a continuous match between <P' and its gradient 
at the bottom of the envelope. This is expressed as 



v = m. 



(25) 



Two boundary conditions must be imposed at the sur- 
face. The first one involves the lagrangian perturbation of 
the pressure dP. As for the radial case, the default choice is 
deduced from the requirement of regularity of the solution 
when P vanishes at the surface. In our variables, this condi- 
tion reads 



0) 



-4 

x^ 



1) q 



1) 



(xP-x 



U 



x^l 



= 0. 



-+x\Z 

X x^ 



(26) 



The user can however choose to impose the more usual con- 
dition 



5^ = or Z = 0. 



(27) 



The second boundary condition ensures the matching of 
and its gradient with the regular solution of the Laplace 
equation outside the star. 



y + (^+i)c/ = o. 



(28) 



Our choice of variables and the way the differential equa- 
tions are written call for three remarks. 
(1) The inclusion of a term in Y in the definition of V (equa- 
tionflSl) ensures that the boundary condition ( |28]) is still valid 
in the (unphysical) case of a non vanishing density at the sur- 
face of the model. 



(2) The use of the lagrangian perturbation of the pression re- 
sults in a better precision in the external layers. 

(3) We do not use the Brunt- Vaisala frequency n in the co- 
efficients of the equations, often badly computed by stellar 
evolution codes. 

However, we use non independent functions pij), m{r) and 
P(r). Troubles may stem from their possible inconsistencies. 
Maybe it would have been wiser to let OSC compute m and 
P from p. 

The solutions computed by OSC are normalized in such 
a way that 



/ 



{cP + l(l+\)b^\dm=MPp 



3.3 Mode classification 

The radial modes owe their existence to the compressibil- 
ity of the stellar material (acoustic or pressure modes). The 
mode with the lowest frequency is called the fundamental 
mode, it has no node in the displacement (except at the ori- 
gin) and is called p\. By order of increasing frequency and 
number of nodes, we have then the first harmonics (one node, 
p-i), the second harmonics (two nodes, pT,), . . . 

For non radial modes, the situation is more complicated. 
For each I, we have a spectrum of 7?-modes. They are of the 
same nature as the radial modes, owing their existence to the 
compressibility of the stellar material. They are numbered 
p\, p2, ...by order of increasing frequency. If the stellar 
model has a radiative zone, it has a spectrum of g+-modes. 
Their frequencies are lower than those of the p-modes and 
have an accumulation point at zero. They are numbered g^, 
g2, ... by order of decreasing frequency. They owe their ex- 
istence to the buoyancy force. For each value of ^ > 1, there 
exists one /-mode, with its frequency between those of the 
p-modes and the g-modes. This mode does not disappear 
when the stellar material is incompressible nor when the 
buoyancy force is zero. When the star has a convective zone, 
another spectrum appears, the g~ -modes. They have an ex- 
ponential temporal behaviour (their frequencies are imagi- 
nary) and are associated with convection. We neglect them 
in the following discussion. 

In OSC, we use an integer to denote the type and or- 
der of a computed mode: n for p„, —n for and for 
/. In the Cowling's approximation (the Eulerian perturba- 
tion of the gravitational potential is neglected), the order of 
a given mode can be easily deduced fro m the behaviour of 
the vertical and horizontal di splace ments dScuflairel Il974bl: 
iGabriel and Scuflairel [l979L Il98(lh . In our case, the mode 
number obtained in the same way is generally correct for 
models which are not too evolved. But as the condensation 
(measured by the ratio pc/p) increases with the age of the 
model, it can just be considered as a clue and finally looses 
any meaning. The algorithm described by iLed (Il985h gives 
a clue to the mode number (also computed by OSC) which 
keeps its utility a bit longer. But when the condensation of 
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the model is really too high, the only reliable identifica- 
tion method consists in the computation of a large number 
of contiguous modes, up to the asymptotic domain, where 
the implemented algorithms continue to give reliable mode 
numbers. 

Though the order of the mode cannot be obtained safely, 
its parity can and is provided by OSC. 



3.4 Influence of rotation 

The rotation of the star removes the degeneracy of the non 
radial oscillation frequencies. If a^^ denotes the frequency 
in the absence of rotation, a slow solid rotation with angular 
velocity Q slightly alters the frequencies in the following 
way 

(ykem = c^u + mPkt^^ , (29) 
with 



{b^ + 2ab)dm 



[a^+e{£+l)b^]din 



(30) 



When the angular velocity depends on the radius, the al- 
tered frequencies may be written 

Okem = + mj Ku {x)Q. [x) dx, (31) 

where the kernel Ku{x), computed by OSC, is given by 

pr^ [a^+ii£+l)b^-2ab-b^] 
KuXr) = f — : — ■ ■ . (32) 



j pr^[a^+l{£+\)b'^] dr 



3.5 Physical description of the modes 

For low order modes, specially for evolved models, the phys- 
ical charac teristics o f a mode is ri ot tightly linked to its g 
or p label (' Scuflairel ll974aL Il980h . OSC outputs different 
indexes allowing the user a quick analysis of the physical 
behaviour of a computed mode (gravity or pressure wave, 
trapped mode, . . . ). 



4 Technique of solution 

4.1 Interpolation of the model 

The grid of points used for the computation of the model 
is rarely appropriate for the computation of oscillations, as 
the eigenfunctions can exhibit rapid spatial oscillations in 
regions where the variables describing the model are well- 
behaved. As the oscillatory behaviour of eigenfunctions is 
easy to foresee, we interpolate the model before any oscil- 
lation computation, increasing the number of points where 
they will prove necessary. The interpolation method we use 
preserves the continuity of the first derivatives. 



4.2 Difference equations 

We have adopted a difference equation scheme of the fourth 
order. That is why we do not need to use Richardson ex- 
trapolation method to increase the precision of the eigenfre- 
quency. 

Our difference scheme rests on the following identity 
satisfied by any vector function y(ji;) with continuous deriva- 
tives up to the fifth order. 



h 



^' 2^' 12^' "^'+'' 



(33) 



where y, and y/_^i are the values of y{x) at points x, and x,+i 
and h = jc,+i — x,. If y is a solution of the linear differential 
system 

dl 
dx 

identity ( [33] | may be written 



■-A{x)y, 



(34) 



{i + ^«, + ^A}y, = { 



= S 1 - 2«'+i 



y<+i 

,5 



with 

a=A, 



dA 
dx 



+0(e), (35) 

(36) 
(37) 



The difference equations are easily obtained from the above 
equations, neglecting the term 0{h^). At the centre, certain 
coefficients of the matrix A are singular and a slightly differ- 
ent treatment is needed. The matrix can be written as 



A{x) = -B{x) 



(38) 



with all the odd order derivatives of matrix B vanishing at 
X = 0. It is clear that the regularity of the solution requires, 
at X = 0, 



By = 0. 



(39) 



It is worth noticing that the rank of B(0) is lower than its 
dimension and that equation (139b gives the right number of 
boundary conditions (1 for the radial case and 2 for the non- 
radial one). The matrices a and j3 assume the following dif- 
ferent forms at the centre (index 0). 



a=0, 



j3=(2-Bo)-' 



d^B\ 
dx^ " 



(40) 
(41) 



4.3 Inverse iteration method 

After the discretization of the differential equations we are 
left with an algebraic eigenvalue problem where the eigen- 
value is /I = (U^. OSC uses the inverse iteration method. It is 
a powerful tool in linear eigenvalue problems. It is described 
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in the b ook of WilkinsonI (I1965L chap 9, sect 47). It has been 
used bv lKeelevl(ll9'77l) to compute stellar radial nonadiabatic 
oscillations. The principle of the method is easy to describe. 
Consider the following eigenvalue problem 

(A-AS)y = 0. (42) 

The method is generally exposed with a unit matrix in place 
of B . Suppose that we know an approximation Ao of an eigen- 
value A. Starting with an arbitrary vector yo, we build the 
sequence 

y„+i =-(A-AoS)-'Sy„ n = 0,l,2,... (43) 
Then, 

Y/i * y«+i 



A - Afl . 



(44) 



yn+l • Yn+l 

The convergence is quite fast. Practically, the y„ must be 
normalized at each step of the computation to avoid over- 
flow. Moreover, these normalized y„ tend to the eigenvector 
associated with A . 

It is true that, in the case of nonradial oscillations, the 
problem to solve is not exactly in the form of equation (1421) . 
But it is put in the right form if we write 

A=Ao + ziA (45) 

and linearize with respect to the correction AX . The solution 
is then obtained in a few iterations of this process. 



5 Use of the program 

OSC is written in Fortran and can be used either as a stand- 
alone program or as a library of subroutines that the user 
calls from a main program of his own. The stand-alone pro- 
gram is in fact just a user interface to the library. It accepts 
instructions from the standard input (or from a command 
file) and prints all kind of information to the standard out- 
put. At the request of the user, the eigenfunctions can be 
saved to a file. For heavy work, the user had better write 
his own Fortran main program and call the routines of the 
library. He gains a better control on the computation and ac- 
ces to results not available otherwise (the rotation kernels for 
a r-dependent angular velocity Q, for instance). 



6 Applications 

OSC is routinely used in our group in Liege and by members 
of the Belgian Asterosesismology Group (BAG) for seis- 
mic studies of solar-lik e pulsators such as, e.g., a Ce n A-i-B 
dlhouLScuflaire et al.', '2003'; 'Miglio and Montalban', '2005^ 
and o f classical /3 Cephei variables ( Aerts et al., 2003; Thoul, 
2003; Dupret et al.LI2004);IAusseloos et al.LI2004) ; IThirion and 



2006.; Briquetet al.r2007h . It is worth mentioning that in 



these studies, we obtained indications on the internal rota- 
tion of the j3 Cephei stars HD 129929 and 6 Ophiuchi. 

The adiabatic frequencies computed by OSC are also 
used as first approximations by the program MAD which 



cor nputes the non adiabatic oscillations of stellar models. 
See lDupretl(l200lh for a description of the code. 

The Liege oscillation code has also taken part in the 

work a nd code comparisons realized within the Corot/ESTA 

group dMontalban and ScuflaireLl2007l : lMontalban. Lebreton et al. 
I2OO7I) . 



7 Discussion 

There is not any general agreement on what the outer bound- 
ary condition on the perturbation of pressure should be and 
different boundary conditions are implemented in the exist- 
ing codes. In our opinion, the precise choice of the outer 
boundary condition does not matter so much, as, in any way, 
the oscillation is far from adiabatic in the very external lay- 
ers of the star. Moreover. , Dzhalilov et al. I (l2000h have shown 
that in the Sun, waves in the frequency range v ~ 2 — 1 mHz 
may reach the chromosphere-corona transition regions by 
means of a tunneling through the atmospheric barrier. When 
comparing with observations, a reasonable strategy is thus to 
use expressions of frequen cies which are insensitive to the 
very external ste llar layers dRoxburgh and Vorontsovl |2003| ; 
lRoxburghl . l2005h . 
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